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Abstract 

In their way to/from turbulence, plane wall-bounded flows dis- 
play an interesting transitional regime where laminar and turbulent 
oblique bands alternate, the origin of which is still mysterious. In 
line with Barkley's recent work about the pipe flow transition involv- 
ing reaction-diffusion concepts, we consider plane Couette flow in the 
same perspective and transform Waleffe's classical four- variable model 
of self-sustaining process into a reaction- diffusion model. We show 
that, upon fulfillment of a condition on the relative diffusivities of its 
variables, the featureless turbulent regime becomes unstable against 
patterning as the result of a Turing instability. A reduced two-variable 
model helps us to delineate the appropriate region of parameter space. 
An intrinsic status is therefore given to the pattern's wavelength for 
the first time. Virtues and limitations of the model are discussed, 
calling for a microscopic support of the phenomenological approach. 

PACS: 47.54.-r (Pattern Formation) ; 47.27.Cn (Turbulent flows. Tran- 
sition to turbulence); 47.20.-k (Flow instabilities. General) 

1 The context 

Patterns are currently observed in continuous media driven out of equilib- 
rium PQ, natural convection being an emblematic case. Pattern formation 
is indeed often an obliged stage in the transition to turbulence. In this re- 
spect, wall-bounded shear flows are systems of great theoretical and practical 
interest. Contrasting with bulk shear flows that become turbulent at low 
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Reynolds number^] in a progressive globally supercritical way, wall-bounded 
flows may remain laminar at pretty high Reynolds numbers under smooth 
enough conditions but experience a direct, discontinuous transition to turbu- 
lence at moderate shear under large enough perturbations, either natural or 
intentional. As a result the transition is subcritical, with hysteresis in some 
range of R called the transitional regime. When confinement effects are weak 
enough, the transition to/from turbulence leaves the realm of chaos theory 
to take on spatiotemporal features. This hysteresis is then responsible for the 
organized coexistence of laminar flow and turbulence in spatially separated 
domains, the patterns of interest here. 

In pipes driven by constant pressure head or mass flux, the transitional 
regime involves chaotic puffs that become turbulent slugs at higher R, see 
e.g. [2] for recent results and references. In plane flows, turbulent spots 
of limited extend [3] can develop to form bands, alternately laminar and 
turbulent. By now, this phenomenon has been observed in several systems, 
Taylor-Couette and plane Couette flow |U El [6] , torsional Couette flow [7] , 
plane Poiseuille flow [8], etc. Being free of global downstream advection 
plane Couette flow, with its pattern at rest in the laboratory frame [BJ, re- 
mains the simplest example. In all uniformly turbulent regime called 
featureless [5J is observed when R is sufficiently large. This occurs beyond 
a threshold value called R t , while below some global stability threshold R g 
laminar flow is always recovered, possibly only at the end of a long turbu- 
lent transient [9]. For plane Couette flow, defining R = Uh/u, with ±U the 
speeds of the counter-translating parallel plates, 2h the gap between them, 
and v the kinematic viscosity of the sheared fluid, one finds R g ~ 325 and 
Rt « 405-415 [SI El QU H]. 

The featureless turbulent regime is well understood in terms of the self- 
sustaining process (SSP) put forward by Waleffe et al. [12] within the Mini- 
mal Flow Unit (MFU) framework [13] • This mechanism was further studied 
by Waleffe in [T4] where it was implemented as a four- dimensional ordinary 
differential system here called Wa97. On the other hand, the emergence of 
bands out of featureless turbulence when R is decreased below Rt has only 
received a phenomenological description in terms of amplitude equations [6J . 
A prior attempt by Hayot and Pomeau [T5] stayed unable to predict any non- 
trivial modulation wavelength while pointing to its possible physical origin. 
To our knowledge, the detailed mechanism however remains unclear and, in 
particular, out of reach of conventional stability analysis of Navier-Stokes 
equations where turbulence is treated by simple closure assumptions [15] . 

1 Comparing shear effects to viscous dissipation over some reievant scale, the Reynolds 
number (R in the following) is the natural control parameter. 
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Recently Barkley introduced an interesting model in which pipe flow is 
considered as an excitable medium described by a reaction-diffusion sys- 
tem [IT], convincingly accounting for the laminar-turbulent dynamics in 
that case [2]. Its extension to plane Couette flow [18] restores the built-in 
upstream/downstream symmetry but involves additional phenomenological 
couplings. Moreover, the excitable character of the local dynamics biases the 
role of laminar and turbulent states in favor of the former, which is reason- 
able only for the laminar-to-turbulent transition close to the lower threshold 
R g . In order to understand the origin of patterning close to upper threshold 
Rt, one would rather want to give the featureless turbulent state a forefront 
role, while staying within the reaction-diffusion framework put forward by 
Barkley. We precisely take this option by considering Wa97 as an appropriate 
starting point. 

In this model, the turbulent state is featured by the upper-branch fixed 
point which is stable when R is large enough. We shall transform the origi- 
nal ordinary differential system into a partial differential system by making 
the degrees of freedom depend also on a space coordinate x, and letting 
them diffuse along that direction. Decreasing R below a well-defined critical 
value Rt (to be identified with the upper threshold Rt), we shall observe 
an instability of the uniform state against a spatial modulation that will be 
interpreted as a standard Turing process [U [19]. This interpretation, giving 
an intrinsic meaning to the pattern's wavelength, will straightforwardly de- 
rive from a reduction of the four-variable model to a two-variable model that 
can be solved by hand. Implications (and limitations) of our findings will be 
discussed next. 

2 The model and some results 

Model Wa97 implement the SSP in the form = J-"(Y; R) where Y is 
a four-component array (M, U, V, W), each variable having a clear physical 
meaning. In a few words, turbulence results from the interplay of the mean 
flow M and streamwise vortices with amplitude V which generate pertur- 
bations called streaks with amplitude U. The streaks are unstable against 
some perturbation W which regenerates the vortices V further distorting M 
via U. Details can be found in [H] to which we refer. 

The very same mechanism operates in pipe flow. The simplified model 
proposed by Barkley [17] involves only two equations for two variables. The 
first variable, it, strictly corresponds to the mean flow M, and the second 
one, q, typifies the turbulence intensity, which can here be identified with 
W. Barkley introduces a coordinate x along the pipe and let variable q 
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diffuse along x through a term d xx q in its own governing equation. The 
upstream/downstream symmetry is broken in the equation for u through a 
term d x u. Global mass advection at some speed U (not to be confused with 
our streak variable) is added. It plays a cosmetic role in the pipe case but 
is crucial to the head-to-tail coupling of two otherwise identical models in 
the Couette case. Such a difficulty is here avoided by letting variables in 
Wa97 diffuse symmetrically along a coordinate x in the direction expected 
for the modulation of turbulence intensity (the pseudo-spanwise coordinate 
z' introduced by Barkley and Tuckerman [10j): 



By assumption the typical scale along x remains unspecified but has to be 
large when compared to the local scales involved in the SSP. System ([ifQ 
will be called 'model Wa97RD' with 'RD' for 'reaction-diffusion'. 

As to the reaction part, the value of coefficients {ay, o~y} with Y = M, U, 
V, or W, explicitly given in [14] . are not interesting in themselves. It however 
warrants to be noted that coefficients accounting for viscous dissipation are 
in the form ay = k y /R where Ky is an effective wavevector associated with 
variable Y, and that nonlinearities in model ([lf|4]) conserve the energy defined 
as E = \{M 2 + U 2 + V 2 + W 2 ) in the same way as the advection term in the 
Navier-Stokes equation preserve the kinetic energy. For what follows, we only 
need to know that the original system has two nontrivial fixed points in 
addition to the trivial fixed point M = 1, U = V = W = which corresponds 
to the linearly stable laminar regime. With the parameters chosen by Waleffe, 
the pair Y^) exists for R > R sn = 104.85, where subscript 'sn' stands for 
'saddle-node'. The so-called lower-branch solution Y* - ) is always unstable, 
while the upper-branch solution Y^ + ) representing the turbulent regime is a 
focus, stable for R > R^ = 138.06 and unstable below ('H' for 'Hopf'). 

Four diffusivities Dy have been introduced, with a priori different values. 
This number can be reduced to three by appropriate rescaling of the space 
coordinate x, i.e. by setting D M = 1, but this still leaves us with three 
independent parameters. 

The stability of the featureless turbulent regime Y(x, t) = Y^ + ) against 
space dependent infinitesimal perturbations is analysed by inserting Y(x, t) = 
Y(+) + Y g exp(st + iqx) in Wa97RD. An instability develops when the real 
part a(q; R) of the complex growth rate s(q; R) = a(q; R) ± iu(q; R) is 
positive. The Turing mechanism works when diffusivities of the 'species' 



d t M + n.y.U 
d t U + a v U 
d t V + a v V 

d t W + a w W 



D M 8 XX M + a M W 2 - o-ulIV + a M 
D v d xx U -a w W 2 + auMV 
D v d xx V + a v W 2 

D w d xx W + a w UW - a M MW - a v VW 



(1) 
(2) 
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(4) 
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in presence, here the velocity components, have sufficiently different mag- 
nitudes [19|. It is characterized by u = and thus generates a pattern 
termed stationary, with some finite wavelength 2Tr/q belonging to an unsta- 
ble band < g m i n < g max - The threshold conditions (Rt, q c ) are defined by 
a(q c ; R T ) = and d q a(q c ; R T ) = 0. 

Before going further, we need an educated guess to fix the value of the 
diffusivities, all other coefficients being given, in order to check whether the 
Turing mechanism gives a possible explanation to the laminar-turbulent band 
formation. Noticing that M and W are the variable closest to those con- 
sidered by Barkley, before considering the general case, we first examine 
what happens when variables U and V are just enslaved to M and W, i.e. 
Djj = Dy = 0, a case which can be solved by hand: 

Anticipating a Turing instability, recalling that u(q; R) = for this mode 
and that perturbations near the threshold are slow, a(q; R) ~ 0, we eliminate 
variables U and V adiabatically by assuming that dtU and dtV are negligible 
when compared to all other terms in (2j3) with Djj = Dy = 0. This yields 
the effective system: 

d t M + a M M = d xx M + a M W 2 + au(Tv(Tw W A - ^MW A + a M ,(5) 

ajjdy otuaty 

(2 2 \ 

— + — ) w 3 
au a v J 

+ °^^MW\ (6) 
auay 



with D := D\y- By construction, system (5|6) has the same fixed points as 



|lj-[4j) and preserve their stability characteristics as long as no intrinsic fre 
quency shows up at threshold. So, the putative Turing mode is not affected 
by the reduction. In contrast, the Hopf instability threshold i?n happens to 
be moved from 138.06 down to R' n = 123.62. The standard approach [19] 
straightforwardly gives the conditions that D must fulfill for a Turing insta- 
bility to develop. Linearization of system (5|6) around fixed point for 



pertubation modes in the form Y q exp(st + iqx) yields an eigenvalue problem: 

(s-g MM + q 2 )M -g MW W = 0, (7) 
-g WM M + (s - g ww + Dq 2 )W = 0, (8) 

where coefficients gyy are easily obtained by explicit computation. The 
characteristic equation of this system is a quadratic polynomial in s with 
coefficients depending on q 2 . The uniform state 'q = 0' must be stable, 
which imposes: 

guu + gww < o , 9mm9ww — gMwgwM > o , (9) 
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Figure 1: As R is decreased at given D, Wa97RD experiences a Turing 
instability for R < i?x given by the curve, with critical wave-vector q c , the 
value of which is indicated for D = 0.05, 0.10, 0.15. Rt = Rn = 138.06 at 
D = 0.1744 with q c = 0.6971. 

When q ^ 0, the sum of the roots £ = guu + Qww ~ (1 + D)q 2 remains 
negative, which forbids any oscillatory instability. A stationary instability 
then develops when the product of the roots 

n = Dq 4 - {Dg MM + g W w)q 2 + 9mm9ww ~ 9mw9wm (10) 

is negative at given D > for some real value of q, hence q 2 > 0, which 
implies Dg MM + gww > 0. The condition that n(g 2 ) = has one or two 
roots reads 

9ww-D 2 + 2(2gMw9wM — guMgww)D + gjVw — > (H) 

which is a condition on D at given gyy- The equality corresponds to the 
double root at threshold, q 2 = q 2 , hence the additional relation dll/d(q 2 ) = 

2Dq 2 = Dg MM + g ww . (12) 

Figure [l] display the result of the threshold conditions {Rt, D} for Wa97RD 
as a curve in the (R, D) plane, the control parameter R being hidden in the 
expressions of the coefficients. The curve itself is obtained from the reduc- 
, which is legitimate since the Turing mode is stationary, but the 



6 




Figure 2: Growth rate a of perturbations to the upper-branch fixed point 
of model Wa97RD as a function of q for D = 0.15 and several values of R 
around the threshold Rt = 143.4. 



unstable domain has been limited to its left by the true condition R > i? H , 
instead of the approximate condition R > R' H stemming from Q. The diffu- 
sivities of the two competing 'species' have thus to be sufficiently different, 
that is D smaller than the critical value computed from condition (11); see 



also §14.3 in [19]. So, for system Wa97RD with D M = 1, D u = D v = 0, and 
Dw '■— D < 0.1744, a Turing instability develops at some threshold > i?n 
and is thus encountered first as R is progressively decreased from large val- 
ues. i?T is seen to increase rapidly as D decreases, being larger than 200 for 
D < 0.05. Figure [2] with D = 0.15, slightly smaller than the limiting value, 
further illustrates the linear stage by displaying the real part of the growth 
rate a(q; R) for different values of R. Two branches appear, the one at small 
q, with u(q; R) ^ not shown, is easily identified as arising from the Hopf 
bifurcation present in the original model. In contrast, the branch at large q 
is stationary, u(q; R) = 0, and corresponds to a standard Turing instability: 
the critical conditions for that example are Rt — 143.4 and q c ~ 0.74. 

A systematic numerical resolution of the full four- dimensional problem 
generalizing system (7|8) has been performed for relative diffusivities Dy/Dm, 
Y = U,V, W, in the form a n , with a = 10 a2 and n G [—15, 15], thus span- 
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Figure 3: Turing unstable domains in parameter space (Dy,Dyi) for dif- 
ferent values of Dy- (a) in the (Dy, Dw) plane for several values of Du] 
(b) in the (Du,D v ) for several values of D w . a = 10 2 1.5849 so that 
a - 15 = 10~ 3 . 

ning the range [0.001, 1000] regularly on a logarithmic scale, D M — 1 fixing 
the scale for coordinate x. The existence of the Turing instability appears 
quite robust, as understood from Figure [3] which displays some significant 
results in the case Dy < 1. It appears that variables M and U on the one 
hand, V and W on the other hand, play on different grounds, and that the 
diffusivities of variables in one group have to be significantly different from 
the diffusivities of the variables in the other group for the Turing mode to be 
relevant. These features are illustrated in the two panels of the figure that 
display isolines Rt = Rn in (Dy,Dy>) planes for the set of Dy considered. 
In this representation, the Turing instability preempts the Hopf instability 
at given Dy when the point corresponding to the values of Dy and Dy of 
interest are in the lower left corner of the panel, below the line labelled by 
Dy"- Panel (a) relative to variables Dy and Dy/ shows that, whatever Du, 
there is no Turing mode if Dy/ > 0.1744 and Dy > 0.0320 and that the larger 
Du the smaller the unstable domain. Clearly, the unstable domain extrapo- 
lates to a well defined region of the (Dy, Dy/) plane, when Du —> 0, which 
is consistent with the results shown in Fig. [T] dedicated to Du = Dy = 0. 
Panel (b) again illustrates the condition Dy < 0.0320 but now shows that 
there can be an instability for Du = 1 (= Dm) if Dy/ < 0.0158. In the fol- 
lowing we shall consider Du = Dm = 1 and Dy = Dy/ := D, in which case 
D < 0.0100 is necessary for the Turing instability to develop. When Du > 1, 
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a qualitatively similar situation holds, which is easier to analyze by scaling 
x using Djj rather than Dm'- the instability is again found when Dy,w are 
small enough and the limit Dm — > behaves like the limit Djj — > in the 
previous case. Quantitative differences however remain because equations 
(II]) and ([2]) are of course not freely exchangable. 

The fate of modulations for R < is also of interest. Here, solutions to 
the full time-dependent nonlinear problem ([lJQ have been obtained by nu- 
merical simulation in a domain of length L = 50, using Neumann boundary 
conditions and a standard finite-difference approach, second order in space 
and time, dealing with the diffusion term implicitly by a Crank-Nicolson 
scheme and the nonlinear interactions explicitly by an Adams-Bashforth 
scheme [20J. Small periodic perturbations around the upper-branch fixed 
point were introduced with a given integer number of cosine arcs, allowing 
us to vary the wavevector per steps Sq = ir/L. The solutions were obtained 
first for R = 140, and next the branches by continuation. Results are pre- 
sented in the form of a bifurcation diagram relating the Reynolds number to 
the amplitude of the steady-state solution defined as the distance to laminar 
flow: A = L- 1 /[(l - Mf + U 2 + V 2 + W 2 } dx. 

In the present work, we mainly consider two extreme cases: the one 
solvable by hand at the linear stage, Du = Dy = 0, and another one with 
Du = Dm large (= 1) and Dy = Dw small (= 0.004). The bifurcation 
diagram corresponding to Du = Dy = and Dw '■— D — 0.15 is shown 
in Figure |4| Besides the amplitude of the non-modulated solutions, the 
figure displays the amplitudes of steady-state solutions with various spatial 
periods as functions of the Reynolds number R^ The bundle of branches 
displayed in the figure corresponds to the whole set of equilibrium solutions 
emerging from initial conditions constructed as described above. Each branch 
is disconnected from the upper-branch base state, and each is terminated by 
two saddle-node bifurcations, one at the high-i? end where the modulated 
state returns to the featureless state, and one at the low-i? end where it 
decays towards the laminar state (A = 0). Remarkably enough, nonlinear 
modulated states can be followed not only below i?n but also well below _R sn , 
when the original model has lost its nontrivial solutions. This is a possibly 
surprising but quite nontrivial effect of the introduction of large scale space 
dependence when passing from Wa97 to Wa97RD. 

In some cases, when the wavevector is large, the upper end point corre- 

2 In fact, these branches could also be reached directly even for R < Rh provided that 
the initial condition be prepared sufficiently close to the upper-branch fixed point. This is 
because the Turing mode has a larger growth rate and develops faster than the Hopf mode 
so that it saturates before the system has a chance to decay via the uniform time-oscillating 
subcritical mode [UIH]. See Fig. [2] 
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Figure 4: Bifurcation diagram of model Wa97RD for D M = 1, D v = D v = 
0, and Dyy = D = 0.15 . Quantity A is given as a function of R for all 
wavevectors q accessible at steady-state in a domain of length L = 50, as 
indicated in the legend. 

sponds to a bifurcation toward a solution with a smaller wavevector, closer 
to the center of the unstable wavevector interval, as expected for an Eckhaus 
instability pQ. In this respect, long wavelength solutions (q small) are much 
more robust than short wavelength ones (q large): we have not been able to 
find equilibrium solutions with q close to the upper bound of the unstable 
wavevector range (Fig. |2|, here for q > 0.8168, while solutions with q close 
to the lower bound, here q = 0.4398, could easily be observed down to very 
low values of R. 

Figure [5^, displays a typical solution obtained at steady-state for D = 0.15 
and R = 135 (< Ru). The laminar state corresponding to M — 1 and 
(U, V, W) = (0,0,0), the laminar-turbulent alternation is easily identified 
with turbulent (laminar) bands associated to the minima (maxima) of M 
and the maxima (minima) of U, V, W. Diffusion implies smooth variations 
of M and W, while the more strongly anharmonic dependence of U and V, 
reaching very low levels inside the laminar regions, is due to their enslaving 
to M and W through nonlinear expressions [T4] . 
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Figure 5: (a) Profile of a stable saturated solution with wavevector q = 
0.6283 at R = 135 for D v = D v = 0, and D w = 0.15. (b) Typical profile of 
the pulse solution at large Djj = Dm and small Dy = Dw for R = 135. 



The interplay between diffusion and nonlinearity is not a trivial mat- 
ter since, when considering our second extreme case, Djj = Dm = 1 and 
Dy = Dw = 0.004, following the same protocol as above, we obtain solu- 
tions that are rapidly driven toward a similar manifold of periodic states, with 
comparable amplitudes, but any of these states is only a long-lived transient 
toward the pulse solution illustrated in Figure [5]d after a cascade of insta- 
bilities progressively reducing the number of oscillation periods. This pulse 
solution is systematically obtained whatever the wavelength of the initial 
condition. It is stable over a very wide range of R: as R is decreased it is ob- 
served down to R ~ 112 below which it decays to the laminar state whereas, 
upon increasing R, it stays localized up to about R = 400 above which it 
expands, triggering the invasion of the upper-branch featureless state. 



3 Discussion 

Up to now, the emergence of laminar-turbulent patterns in transitional wall- 
bounded flows has not receive any clear-cut explanation [TB] . Elaborating on 
ideas put forward by Barkley [TTl HE] who interprets the laminar-to-turbulent 
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transition in physical space as the result of excitatory-refractory behavior 
common in reaction-diffusion processes [19], we have shown that a model 
introduced by Waleffe [14J to describe the local sustainment of turbulence, 
once appropriately converted into a reaction-diffusion system, could account 
for the development of a pattern similar to what is observed experimentally 
in plane Couette flow at the turbulent-to-laminar transition [6j. This under- 
standing in terms of a Turing instability of the featureless turbulent regime 
(Fig. [l] & [2]) points to a possible generic origin of the phenomenon, and ac- 
cordingly could also apply in other similar situations [H [5] : Provided that 
large scale perturbations evolve on sufficiently different spatiotemporal scales, 
as a result of effective diffusivities of sufficiently different magnitudes, we have 
shown that infinitesimal modulations of the turbulence intensity around the 
featureless state, here represented by the upper-branch fixed point of Wal- 
effe's model, ends in a nontrivial patterning in the transitional range. 

Let us first note that the transformation of a local model of SSP into a 
reaction-diffusion system via the phenomenological introduction of effective 
diffusion terms is not as arbitrary as it might seem since it aims at accounting 
for large scale modulations of the SSP intensity in much the same way as the 
eddy viscosity helps us at managing the effects of turbulent fluctuations on 
the mean flow in standard turbulence theory. 

Next, several features retrieved from experiments are satisfactorily ren- 
dered within the model as it stands. The grouping of variables relevant to the 
SSP mechanism in two distinct sets, {M, U} and {V, W}, was not obvious in 
advance, nor the ordering of the diffusivities (Fig. [3]), though this property 
is in line with longer coherence lengths for the mean flow M and streaks U, 
than for the streamwise vortices V and the streak instability mode W, as 
expected from observations (22]. 

At the nonlinear stage, results are less satisfactory since the instability 
was found discontinuous at R? while its seems to be continuous in experi- 
ments [B]. The fact that, in all cases, nonlinear solutions can be found at 
values of R where the featureless state is unstable against uniform modes 
(R < Ru) and, in some cases, well below R sn (Fig. El), is however a strong 
indication that large-scale spatiotemporal couplings profoundly modify the 
small-scale picture gained using the MFU assumption [T3], e.g. the search 
for invariant temporal solutions within the framework of dynamical systems 
theory [23J. 

Despite the appealing features of the Turing instability concept, and es- 
pecially the internal nature of the patterning mechanism, vis. Eq. (12) fixing 
the critical wavelength, the weak point of our approach lies in its phenomeno- 
logical nature, implying an indeterminate direction for the space coordinate 
and an arbitrary scale along that direction. In view of the prediction of quan- 



12 



titative properties of the pattern, orientation and wavelength, some support 
of the reaction-diffusion formalism has thus to be gained from a more ab 
initio (microscopic) approach. An extension of previous Galerkin model- 
ing [21] combined to filtering able to separate large and small scales and its 
adaptation to other cases of great interest for applications [5j [?l E] is under 
development in this purpose. 
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